Temperature dependence of the paramagnetic spin susceptibility of doped graphene 
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In this work, we present a semi-analytical expression for the temperature dependence of a spin- 
resolved dynamical density-density response function of massless Dirac fermions within the Random 
Phase Approximation. This result is crucial in order to describe thermodynamic properties of the 
interacting systems. In particular, we use it to make quantitative predictions for the paramagnetic 
spin susceptibility of doped graphene sheets. We find that, at low temperatures, the spin susceptibil- 
ity behaves like T~ 2 which is completely different from the temperature dependence of the magnetic 
susceptibility in undoped graphene sheets. 
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I. INTRODUCTION 

Graphene is a newly realized two-dimensional (2D) 
electron system that has attracted a great deal of in- 
terest because of the new physics which it exhibits and 
because of its potential as a new material for electronic 
technology 1,2 . It exhibits a large number of new and 
exotic optical and electronic effects that have not been 
observed in other materials 3 . 

When non-relativistic Coulombic electron-electron 
interactions are added to the kinetic Hamiltonian, 
graphene represents a new type of many-electron prob- 
lem, distinct from both an ordinary 2D electron gas (EG) 
and from quantum electrodynamics. The Dirac-like wave 
equation and the chirality of its eigenstates lead indeed to 
both unusual electron-electron interaction effects 4 - - — and 
an unusual response to external potential s 10 ' 11 . 

Spin transport of fermions is central to many fields of 
physics. Electron transport runs modern technology and 
electron spin is being explored as a new carrier of in- 
formation 12 . There has been recent interest in the tem- 
perature dependence of various Fermi liquid properties 13 . 
Technical advances now allow one to measure the temper- 
ature dependence of the thermodynamic and transport 
parameters such as conductivity at finite temperature 
in unsuspende d 14 ' 15 , and suspended graphene sheet s 16 ' 17 
and spin susceptibility in two-dimensional systems. 

The paramagnetic spin susceptibility shows behavior 
similar to the charge compressibility 6 which decreases as 
the interaction increases at zero temperature. This is re- 
lated to the fact that the inverse of the spin susceptibility 
is proportional to the renormalized Fermi velocity. 

Vafeki^ has recently shown that the specific heat of 
undoped graphene sheets presents an anomalous low- 
temperature behavior displaying a logarithmic suppres- 
sion with respect to its noninteracting counterpart as 
~ T/ln(T). Meanwhile, Shcery and Schmalian studied 
both the specific heat and the orbital diamagnetic suscep- 
tibility of undoped graphene by using a renormalization 
group approach™. They stated that the dependence of 
the diamagnetic susceptibility of undoped graphene on 
temperature is quite different from the 2D EGs and it 



behaves like T/|ln(T)| 2 . 

On the other hand, it has been demonstrated in 
RefspiiS (see also Ref3) that doped graphene sheets are 
normal (pseudochiral) Fermi liquids, with Landau pa- 
rameters that possess a behavior quite distinct from 
that of conventional 2D EGs. In addition, it was found 
that^S, at low temperatures, the specific heat has the 
usual normal-Fermi-liquid linear-in-temperature behav- 
ior, with a slope that is solely controlled by the renor- 
malized quasiparticle velocity in doped graphene. 

The temperature correction to the spin susceptibility 
for a 2D EG interacting via a long-range Coulomb in- 
teraction has attracted a lot of interest over a long pe- 
riod of time 21 -. It has been shown that the dynamic 
Kohn anomaly in the density-density response function 
at 2fcp and re-scattering of pairs of quasiparticles lead 
to linear-in-temperature correction to the spin suscepti- 
bility 2 -. Since the static non-interacting density-density 
response function of doped graphene is a smooth function 
at 2fcp and behaves differently from what one has in stan- 
dard 2D EG systems, a linear-in-temperature correction 
to the spin susceptibility does not occur. 

In this work we calculate the temperature depen- 
dence of a spin-resolved dynamical density-density re- 
sponse function of massless Dirac fermions within the 
Random Phase Approximation (RPA) and subsequently 
the Helmholtz free energy -F(T) of doped graphene sheets 
where the chemical potential is non zero. This allows us 
to access important thermodynamic quantities, such as 
the spin susceptibility which can be calculated by taking 
appropriate derivatives of the free energy. We show that, 
at low temperatures, the spin compressibility of doped 
graphene, in contrary to the diamagnetic spin susceptibil- 
ity^, behaves as the inverse square of temperature, solely 
controlled by both the ultraviolet cut-off and graphene's 
fine-structure constant. 

In Sec. II we introduce the formalism that will be used 
in calculating (paramagnetic) spin susceptibility which 
includes the many-body effects in the RPA. In Sec. Ill we 
present our analytical and numerical results for the free 
energy and spin susceptibility in doped graphene sheets. 
Sec. V contains discussions and conclusions. 
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II. METHOD AND THEORY 

The agent responsible for many of the interesting elec- 
tronic properties of graphene sheets is the non-Bravais 
honeycomb-lattice arrangement of carbon atoms, which 
leads to a gapless semiconductor with valence and con- 
duction 7r-bands. States near the Fermi energy of 
a graphene sheet are described by a spin-independent 
massless Dirac Hamiltonian 2 ^ 



H = Vp(T ■ p , 



(1) 



where vp is the Fermi velocity, which is density- 
independent and roughly three-hundred times smaller 
that the velocity of light in vacuum, and er = (a x ,a y ) is a 
vector constructed with two Pauli matrices {a\ i — x,y}, 
which operate on pseudospin (sublattice) degrees of free- 
dom. Note that the eigenstates of H have a definite chi- 
rality rather than a definite pseudospin, i.e. they have 
a definite projection of the honeycomb-sublattice pseu- 
dospin onto the momentum p. 

Electrons in graphene do not move around as inde- 
pendent particles. Rather, their motions are correlated 
due to pairwise Coulomb interactions. The interaction 
potential is sensitive to the dielectric media surrounding 
the graphene sheet. The Fourier transform of the real 
space potential is given by v q = 2ne 2 /eq where e is the 
average dielectric constant between the medium and a 
dielectric constant of the substrate. 

Within this low energy description, the properties of 
doped graphene sheets depend on the dimensionless cou- 
pling constant or graphene's fine-structure constant 



(2) 



ehvp 



As it is clearly seen from Eq. (JTJ), the spectrum is un- 
bounded from below and it implies that the Hamiltonian 
has to be accompanied by an ultraviolet cut-off which is 
defined k c and it should be assigned a value correspond- 
ing to the wavevector range over which the continuum 
model Eq. (|TJ) describes graphene. For later purposes we 
define a gr = 2g v a cc where g v = 2 accounts for valley 
degeneracy, kp = kp(l + (cr) 1 ^ 2 is the spin dependent 
Fermi momentum and kp = (7m) 1 / 2 is the Fermi wave 
number and ep = hv-pk-p is being the Fermi energy with 
n = n 1 ' + nr the total electron density, £ = (rzT — n^) jn 
is the spin polarization parameter (0 < £ < 1 if we as- 
sume that, e. g., electrons with real spin s =t to be 
majority). For definiteness we take fc c to be such that 
7rfc 2 = 2(27r) 2 /.4o, where Ao — 3V3a 2 ,/2 is the area of the 
unit cell in the honeycomb lattice, with ao — 1.42 A the 
Carbon-Carbon distance. With this choice, the energy 
Hvk r = 7 eV and 



2gv 
nAo 



(3) 



The continuum model is useful when k c ^$> kp, i.e. when 
A > 1. Note that, for instance, electron densities n = 



0.36 x 10 12 and 0.36 x 10 14 cm" 2 correspond to A = 100 
and 10, respectively. 

The free energy T = Ta + ^"int , a thermodynamic po- 
tential at a constant temperature and volume, is usually 
decomposed into the sum of a noninteracting term Tq 
and an interaction contribution T^t- To evaluate the in- 
teraction contribution to the Helmholtz free energy we 
follow a familiar strateg y 24 ] 25 by combining a coupling 
constant integration expression for J-- m t valid for uniform 
continuum models (h = 1 from now on), 



■Fint(T) 



N 
~2 



dX 



d 2 q 

(27T) 



2 U <1 



(4) 



with a fluctuation-dissipation-theorem (FDT) expres- 
sion^ 5 - for the static structure factor, 



1 



dtu coth(/3uj/2)Qmx ( p X J(q,uj,T) 



S W (q,T) 

(5) 

Here (3 = (k&T) 1 . Quite generally, two-particle correla- 
tion functions can be written in terms of single-particle 
Greens functions and vertex parts. The RPA approxima- 
tion for J^nt then follows from the RPA approximation 

for Xpp{Qi u ) : 



1 - >>vq(xl{q, t) + xi(q, w, T)) 



(6) 



where Xoil^^) ^ s ^ e noninteracting spin resolved 
density-density response-function in er channel. A cen- 
tral quantity in the many-body techniques is the nonin- 
teracting spin resolved polarizability function Xo (<7j w j T) 
. The problem of linear density response is set up by con- 
sidering a fluid described by the Hamiltonian, which is 
subject to an external potential. The external potential 
must be sufficiently weak for low-order perturbation the- 
ory to suffice. The induced density change has a linear 
relation to the external potential through the noninter- 
acting dynamical polarizability function. This function 
in er channel reads as 

a, ^ v ST f d2k l + ss'cos(6»fc !fe+ q) 
Xo(s, U ,T) = 5% hm ^ J „ 



, (2tt) 2 

s,s>=±" V ' 

n F (e kiS ) - n F (e k + q ,s') 
to + Sk,s - £k+ q ,s' + if] 



(7) 



Here Sk , s — svpk are the Dirac band energies and 
n p{ £ ) = { ex P[/^( e — A*o)] + i s t ne usua l Fermi-Dirac 
distribution function, /1q = /ifi (T) being the noninteract- 
ing chemical potential. As usual, this is determined by 
the normalization condition 



de v{e)n F {e) 



(8) 



where v(e) = g v e/(2irvp) is the noninteracting den- 
sity of states. For T -> one finds ^o( T ) = £ f{ 1 ~ 
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7r 2 (T/T F ) 2 /6(l + crC)}, where T F = £ F /£b is the Fermi 
temperature. The factor in the first line of Eq. ([7]), which 
depends on the angle 6k,k+q between k and k + q, de- 
scribes the dependence of Coulomb scattering on the rel- 
ative chirality ss' of the interacting electrons. 

After some straightforward algebraic manipulations^ 
we arrive at the following expressions for the imaginary, 
3m Xo ) ano - the real, $te Xo> parts of the noninteracting 
density-density response function for uj > 0: 



Xo(?,w,T) = ^ \ B(w F <7- uj)q 2 f(v F q,uj) 
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G^' c7) (q,w,T)-G (Q ' c7) (g,w,r) 
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^8 a ,- + H<*°\q, U ,T) 



(9) 



and 
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FIG. 1: ( Color online) Upper : the real part of the dynam- 
ical response function Ke Xoil^iT) [in units of v{ef)] as a 
0(u - v F q)q 2 f{uj,v F q) Ig*?'^ (q,u>,T) - G { +' a) (q, u, T) Wet ion of q/k F for u) = 2e P , T = 0.1T F and three values of 

-0 < C < 1- Bottom: same as in the upper panel but for the 
imaginary part. 



+ e(v F q~uj)q 2 f(v F q,uj) ~8 a ,- + H^"\q, u, T) 



Here 



ution to the free energy per particle fmt{T) is given by 

I f d 2 q f 1 '■ +o ° 



/int(T) 



f(x,y) 



2 V / ^ 3 " 



(11) 



arctan 



1 

2ri 



2 J (2tt) 2 
Vql^m Xq + Sto Xq] 
l-v g [$le xl + Ke xb) 
d 2 q f 1 dX 



dui coth (Puj/2)x 



and 



exp 



2fc^T 



(2^) 2 y a 

d[^e xl + Me xt] 



coth(/3w pl /2)[3?e xl + ^e xt] 



dui 



(14) 



+ 1 
(12) 



v F qu ± lo\ — 2a^Q 



(13) 

The coupling constant integration in Eq. Q can be 
carried out partly analytically due to the simple RPA 
expression Eq. ([5]). We find that the interaction contri- 



In this equation the first term comes from the smooth 
electron- hole contribution to 3to Xpp 1 while the sec- 
ond term comes from the plasmon contribution; ui^ = 
ui p \(q, T, A) is the plasmon dispersion relation at coupling 
constant A which can be found numerically by solving the 
equation 1 — Af q 5Re Xo (<7; w ; ^) = 0- Note that in a stan- 
dard 2D EG the exchange energy starts to matter little 
for T of order T F because all the occupation numbers 
are small and the Pauli exclusion principle matters little. 
In the graphene case however exchange interactions with 
the negative energy sea remain important as long as T is 
small compared to v F k c /k^ — T F A. 
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FIG. 2: (Color online) Upper: <5/i n t in units of ep as a function 
of spin polarization parameter, £ for a gr = 1 and T = O.ITf. 
Bottom : Numerical calculated SxJ 1 ^) = X^ 1 ( T )-X^ 1 ( r = 
0) in units of £F/nfi% as a function of temperature for a gr = 1 
in comparison with the low temperature approximated ex- 
pression given by Eq. (ll9p . These numerical results confirm 
the validity of our analytic result for <5/mt- In the inset, the 
inverse spin susceptibility scaled by non-interacting one as 
a function of electron density in units of 10 1 
temperature for different a gr values. 



The free energy calculated according to Eq. (fl~4"|) is di- 
vergent since it includes the interaction energy of the 
model's infinite sea of negative energy particles. Follow- 
ing Vafeki^, we choose the free energy at T = 0, /(T = 
0), as our "reference" free energy, and thus introduce the 
regularized quantity Sf = f{T) — f(T = 0). This again 
can be decomposed into the sum of a noninteracting con- 
tribution, Sf (T 0) = - 3 „£ F 7r 2 (T/T F ) 2 Z(C,l/2)/12, 
where Z((, m) = (1 + C) m + (1 — C) m an d an interaction- 
induced contribution Sf int {T) = f hlt (T) ~ f- m t(T = 0), 
which can be calculated from Eq. (T14")) . Note that we 
have f (T = 0) = g v e F Z{(, 3/2)/6. 

The low-temperature behavior of the interaction con- 
tribution to the free energy can be extracted analytically 
with some patience. After some lengthy but straightfor- 
ward algebra we find, to leading order in A, 



FIG. 3: (Color online) Upper: the spin susceptibility (in units 
of the non-interacting spin susceptibility xo) as a function of 
coupling constant for two values of ultraviolet cut-off, A = 10 
and 100. The spin susceptibility decreases with increasing A 
or the coupling constant. Bottom: the same as upper panel 
for three values of temperature for A = 100. 



<5/int(T^0) = £F - 
xZ(C,l/2)lnA + R. T 



a gr [l - a gr £(a gr ) 



(15) 



where the function £(x), defined as in Eq. (14) of Ref. 0, 
is given by £(x) = 128/(tt 2 x 3 ) - 32/(tt 2 x 2 ) + 1/x - 
h(irx/8), with 



h(x) = < 



2x 3 VT 



: arctan 



Vl^ 2 



4cc 3 Va 



:ln 



+ Va?2 



1 



- Vx^l 



for x < 1 



for x > 1 



(16) 

The symbol R. T. in the l.f.s. of Eq. (fT5j) indicates regular 
terms, i.e. terms that, by definition, are finite in the limit 
A — > oo. Eq. (|15p represents the second important result 
of this work. 
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We thus see that Sf int (T -> 0) cx T 2 in Eq. JTSJ im- 
plies a conventional Fermi-liquid behavior with a linear- 
in-T specific heat. We are thus led to conclude, in full 
agreement with the zero-temperature calculations of the 
quasiparticle energy and lifetime performed in Refsi 7 -^, 
that doped graphene sheets are normal Fermi liquids. 

It would be worthwhile obtaining the high-temperature 
dependence of Sf. Since we are always measuring 
energies in units of the Fermi energy, therefore our 
high energy results are relevant to Dirac point physics. 
The undoped limit for us is the limit of vanishing 
the Fermi energy. To obtain the results for undoped 
graphene, let's consider the paramagnetic case. By re- 
placing the Fermi energy with k-gT and therefore, fcp 
with ksT/hv-p in Eq. QjjJ the correction to the free- 
energy is given by 5f- mt (T > T F ,n ~ 0) - T 3 a gr [l - 
a gr £(a gr )]hx(k c /T)/(8g v ). Importantly, this expression, 
apart from a constant^, is coincide with that result ob- 
tained in Eq. 13 by Vafek^. Therefore, we expect that in 
the limit that T 3> Tp and for every £ value, the temper- 
ature dependence of the free-energy correction behaves 
like T 3 ln(fc c /T). 

The spin susceptibility, on the other hand, can be cal- 
culated from the second derivative of the Helmholtz free 
energy and it reads 

i _ i 9 2 [/ (T,c) + /int(r,C)] . ri7 . 

xMi lc =° ' (17) 

where [Ib is the Bohr magneton 27 . 
It is easy to calculate the non-interaction spin suscep- 
tibility and it turns out that 

At low temperature, by using 5/i n t to leading order in 
A, the temperature dependence of the correction to the 
spin susceptibility is thus given by 

^CO = x; 1 (t)- x ; 1 (t = o) 

£ ftt 2 / T 2 \g v 1 

= % } iy- 7 ^ (19) 

where r\ — a gr (l — a gr £,(a gr ))/12g v . It is obvious from 
the expression that Xs{T) cx T~ 2 at low temperature 
limit. This expression represents our important result in 
this work. 



III. NUMERICAL RESULTS 

In this section, we present the most important results 
of the spin susceptibility in doped graphene sheets by 
using mentioned formalism. 

The semi-analytical expressions for 5ie \o (<7i w > T) and 
Qm Xo(9 5 w,T) constitute the first important result of 
this work. In Fig. 1 we have plotted the major part of 



the dynamic response function as a function of q/k-p for 
different values of £. Sharp cutoffs in the imaginary part 
of Xg(g, td,T) are related to the rapid swing in the real 
part of xJ(g,u;,T). These behaviors are in result of the 
fact that the real and imaginary parts of the polarization 
function are related through the Kramers-Kronig rela- 
tions. Importantly, the sign change of the real part from 
negative to positive shows a sweep across the electron- 
hole continuum. It is important to note that there is a 
non-monotonic behavior of Xq(q, 0,T) as a temperature 
dependent originates from a competition between intra- 
and inter-band contributions to this quantity 2 ^ How- 
ever, the spin polarization parameter dependence of the 
Lindhard function is a monotonic behavior at any fre- 
quency. 

Fig. 2 (upper) shows the interaction contributions of 
the free energy as a function of spin polarization for 
T = O.lTp. Our numerical results show that the contri- 
bution of electron-electron interactions in the free energy 
decrease by increasing the spin polarization. It changes 
slightly at low £ values and sharply decreases near the 
ferromagnetic case. The reason is that the exchange and 
correlation energies mostly change around the fully fer- 
romagnetic point. Moreover, the slope of exchange and 
correlation energies with respect to £ around £ = 1 have 
opposite signs 28 and the interaction contribution tends to 
the correlation sign at certain value of £■ Fig. 2 (bottom) 
shows the numerical calculated x i T 1 (r) — X^ X (T — 0) as 
a function of temperature in comparison with that re- 
sult obtained at low temperature and leading order of A. 
We can easily see that those results are very close at low 
temperature and confirm that the spin susceptibility be- 
haves as T~ 2 in this region. In addition, this comparison 
allows us to use the approximated analytical expression 
given by Eq. (|19[) for the temperature correction of the 
spin susceptibility till T < 0.3Tp. Furthermore, we can 
see that the spin susceptibility sharply decreases with in- 
creasing temperature. On the other hand, temperature 
dependence of the spin susceptibility is in contrast to that 
result obtained for the diamagnetic undoped graphene 
sheet. To seek comprehensive study, we have numeri- 
cally calculated Xa/Xs as a function of electron density 
in units of 10 12 cm -2 at zero temperature, T = 0. Our 
results are shown in the inset of Fig. 2 (bottom) and 
show that the spin susceptibility increases by increasing 
the electron density^. 

Finally, we show the spin susceptibility scaled by its 
non-interacting value as a function of the coupling con- 
stant for (upper panel) two values of the ultraviolet cut- 
off and (bottom panel) different values of temperatures 
in Fig. 3. These results are obtained numerically by tak- 
ing the full terms of Eq. (TT^l) . We can clearly see that 
the spin susceptibility increases by increasing the electron 
density ( or decreasing the A values) while it decreases by 
increasing the interactions at certain temperature value. 
The reason is that the exchange contribution term makes 
a positive contribution to Eq. (|18[) . thus tending to re- 
duce the spin susceptibility (with respect to its nonin- 
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teracting value), again in contrast to what happens in 
the standard 2D EG where exchange enhances the spin 
susceptibility. The correlation term instead makes a neg- 
ative contribution to Eq. (j!8[) . thus tending to enhance 
the spin susceptibility. In the 2D EG, correlations tend 
to reduce the spin susceptibility. 

IV. CONCLUSIONS 

We have presented semi-analytical expressions for the 
real and the imaginary parts of the resolved spin de- 
pendence of density-density linear-response function of 
noninteracting massless Dirac fermions at finite temper- 
ature. These results are very useful in order to study 
finite-temperature screening within the Random Phase 
Approximation. For example they can be used to cal- 
culate the spin dependence of the conductivity at finite 
temperature within Boltzmann transport theory . 

The Lindhard function at finite temperature is also 
extremely useful to calculate finite-temperature equilib- 
rium properties of interacting massless Dirac fermions, 



such as the specific heat and the compressibility. For ex- 
ample, in this work we have been able to show that, at 
low temperatures, the paramagnetic spin susceptibility of 
interacting massless Dirac fermions behaves like T~ 2 at 
low temperature. Even though the charge and spin sus- 
ceptibilities behave similarly at zero temperature, their 
temperature dependencies are totally different. We have 
obtained an analytical expression for the spin susceptibil- 
ity in the leading order of cut-off and showed that one can 
use that in the low temperature range for experimental 
access. 

We remark that in a very small density region, the 
system is highly correlated and a model going beyond 
the RPA is necessary to account for increasing correlation 
effects at low density. 
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